library(quantreg)
library(affy)

setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/1")
data_1<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/2")
data_2<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/3")
data_3<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/4")
data_4<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/5")
data_5<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/6")
data_6<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/7")
data_7<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/8")
data_8<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/9")
data_9<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/10")
data_10<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/11")
data_11<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/12")
data_12<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/13")
data_13<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/14")
data_14<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/15")
data_15<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/16")
data_16<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/17")
data_17<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/18")
data_18<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/19")
data_19<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/20")
data_20<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/21")
data_21<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/22")
data_22<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/23")
data_23<-read.table("normalize.txt",sep="\t",head=T)
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/24")
data_24<-read.table("normalize.txt",sep="\t",head=T)

Cen_1F<-apply(data_1,1,mean)
Cen_2F<-apply(data_2,1,mean)
Cen_3F<-apply(data_3,1,mean)
Cen_4F<-apply(data_4,1,mean)
Cen_5F<-apply(data_5,1,mean)
Cen_6F<-apply(data_6,1,mean)
Cen_7F<-apply(data_7,1,mean)
Cen_8F<-apply(data_8,1,mean)
Cen_9F<-apply(data_9,1,mean)
Cen_10F<-apply(data_10,1,mean)
Cen_11F<-apply(data_11,1,mean)
Cen_12F<-apply(data_12,1,mean)
Cen_13F<-apply(data_13,1,mean)
Cen_14F<-apply(data_14,1,mean)
Cen_15F<-apply(data_15,1,mean)
Cen_16F<-apply(data_16,1,mean)
Cen_17F<-apply(data_17,1,mean)
Cen_18F<-apply(data_18,1,mean)
Cen_19F<-apply(data_19,1,mean)
Cen_20F<-apply(data_20,1,mean)
Cen_21F<-apply(data_21,1,mean)
Cen_22F<-apply(data_22,1,mean)
Cen_23F<-apply(data_23,1,mean)
Cen_24F<-apply(data_24,1,mean)

CenF<-(Cen_1F+Cen_2F+Cen_3F+Cen_4F+Cen_5F+Cen_6F+Cen_7F+Cen_8F+Cen_9F+Cen_10F+Cen_11F+Cen_12F+Cen_13F+Cen_14F+Cen_15F+Cen_16F+Cen_17F+Cen_18F+Cen_19F+Cen_20F+Cen_21F+Cen_22F+Cen_23F+Cen_24F)/24
betweenF <- 5*((Cen_1F-CenF)^2+(Cen_2F-CenF)^2+(Cen_3F-CenF)^2+(Cen_4F-CenF)^2+(Cen_5F-CenF)^2+(Cen_6F-CenF)^2+(Cen_7F-CenF)^2+(Cen_8F-CenF)^2+(Cen_9F-CenF)^2+(Cen_10F-CenF)^2+(Cen_11F-CenF)^2+(Cen_12F-CenF)^2+(Cen_13F-CenF)^2+(Cen_14F-CenF)^2+(Cen_15F-CenF)^2+(Cen_16F-CenF)^2+(Cen_17F-CenF)^2+(Cen_18F-CenF)^2+(Cen_19F-CenF)^2+(Cen_20F-CenF)^2+(Cen_21F-CenF)^2+(Cen_22F-CenF)^2+(Cen_23F-CenF)^2+(Cen_24F-CenF)^2 )/23
dataF<-as.matrix(cbind(data_1[,1:5],data_2[,1:5],data_3[,1:5],data_4[,1:5],data_5[,1:5],data_6[,1:5],data_7[,1:5],data_8[,1:5],data_9[,1:5],data_10[,1:5],data_11[,1:5],data_12[,1:5],data_13[,1:5],data_14[,1:5],data_15[,1:5],data_16[,1:5],data_17[,1:5],data_18[,1:5],data_19[,1:5],data_20[,1:5],data_21[,1:5],data_22[,1:5],data_23[,1:5],data_24[,1:5]))

acrossF<-CenF
for(i in 1:nrow(dataF)){
acrossF[i]<-0
}
for(i in 1:ncol(dataF)){
	if(i<6)
	{
		temp<-((dataF[,i]-Cen_1F)^2)/96
		acrossF<-acrossF + temp
	}
	if(i<11 &&i>5)
	{
		temp<-((dataF[,i]-Cen_2F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<16 && i>10)
	{
		temp<-((dataF[,i]-Cen_3F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<21 && i>15)
	{
		temp<-((dataF[,i]-Cen_4F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<26 && i>20)
	{
		temp<-((dataF[,i]-Cen_5F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<31 && i>25)
	{
		temp<-((dataF[,i]-Cen_6F)^2)/96
		acrossF<-acrossF + temp
	
	}
	if(i<36 && i>30)
	{
		temp<-((dataF[,i]-Cen_7F)^2)/96
		acrossF<-acrossF + temp
	}
	if(i<41 &&i>35)
	{
		temp<-((dataF[,i]-Cen_8F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<46 && i>40)
	{
		temp<-((dataF[,i]-Cen_9F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<51 && i>45)
	{
		temp<-((dataF[,i]-Cen_10F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<56 && i>50)
	{
		temp<-((dataF[,i]-Cen_11F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<61 && i>55)
	{
		temp<-((dataF[,i]-Cen_12F)^2)/96
		acrossF<-acrossF + temp
	
	}
	if(i<66 && i>60)
	{
		temp<-((dataF[,i]-Cen_13F)^2)/96
		acrossF<-acrossF + temp
	}
	if(i<71 &&i>65)
	{
		temp<-((dataF[,i]-Cen_14F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<76 && i>70)
	{
		temp<-((dataF[,i]-Cen_15F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<81 && i>75)
	{
		temp<-((dataF[,i]-Cen_16F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<86 && i>80)
	{
		temp<-((dataF[,i]-Cen_17F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<91 && i>85)
	{
		temp<-((dataF[,i]-Cen_18F)^2)/96
		acrossF<-acrossF + temp
	
	}
	if(i<96 && i>90)
	{
		temp<-((dataF[,i]-Cen_19F)^2)/96
		acrossF<-acrossF + temp
	}
	if(i<101 &&i>95)
	{
		temp<-((dataF[,i]-Cen_20F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<106 && i>100)
	{
		temp<-((dataF[,i]-Cen_21F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<111 && i>105)
	{
		temp<-((dataF[,i]-Cen_22F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<116 && i>110)
	{
		temp<-((dataF[,i]-Cen_23F)^2)/96
		acrossF<-acrossF + temp

	}
	if(i<121 && i>115)
	{
		temp<-((dataF[,i]-Cen_24F)^2)/96
		acrossF<-acrossF + temp
	
	}

}
F<-betweenF/acrossF
setwd("A:/data_HTHGU/maqc/affy_plus2_taqman/residualplot/plot4_inputs")
data<-read.table("rma_normalizaed.txt",sep="\t",head=T)
Cen_1<-apply(data[,1:5],1,mean)
Cen_2<-apply(data[,6:10],1,mean)
Cen_3<-apply(data[,11:15],1,mean)
Cen_4<-apply(data[,16:20],1,mean)
Cen_5<-apply(data[,21:25],1,mean)
Cen_6<-apply(data[,26:30],1,mean)
Cen_7<-apply(data[,31:35],1,mean)
Cen_8<-apply(data[,36:40],1,mean)
Cen_9<-apply(data[,41:45],1,mean)
Cen_10<-apply(data[,46:50],1,mean)
Cen_11<-apply(data[,51:55],1,mean)
Cen_12<-apply(data[,56:60],1,mean)
Cen_13<-apply(data[,61:65],1,mean)
Cen_14<-apply(data[,66:70],1,mean)
Cen_15<-apply(data[,71:75],1,mean)
Cen_16<-apply(data[,76:80],1,mean)
Cen_17<-apply(data[,81:85],1,mean)
Cen_18<-apply(data[,86:90],1,mean)
Cen_19<-apply(data[,91:95],1,mean)
Cen_20<-apply(data[,96:100],1,mean)
Cen_21<-apply(data[,101:105],1,mean)
Cen_22<-apply(data[,106:110],1,mean)
Cen_23<-apply(data[,111:115],1,mean)
Cen_24<-apply(data[,116:120],1,mean)
Cen<-(Cen_1+Cen_2+Cen_3+Cen_4+Cen_5+Cen_6+Cen_7+Cen_8+Cen_9+Cen_10+Cen_11+Cen_12+Cen_13+Cen_14+Cen_15+Cen_16+Cen_17+Cen_18+Cen_19+Cen_20+Cen_21+Cen_22+Cen_23+Cen_24)/24
between <- 5*((Cen_1-Cen)^2+(Cen_2-Cen)^2+(Cen_3-Cen)^2+(Cen_4-Cen)^2+(Cen_5-Cen)^2+(Cen_6-Cen)^2+(Cen_7-Cen)^2+(Cen_8-Cen)^2+(Cen_9-Cen)^2+(Cen_10-Cen)^2+(Cen_11-Cen)^2+(Cen_12-Cen)^2+(Cen_13-Cen)^2+(Cen_14-Cen)^2+(Cen_15-Cen)^2+(Cen_16-Cen)^2+(Cen_17-Cen)^2+(Cen_18-Cen)^2+(Cen_19-Cen)^2+(Cen_20-Cen)^2+(Cen_21-Cen)^2+(Cen_22-Cen)^2+(Cen_23-Cen)^2+(Cen_24-Cen)^2 )/23
across<-Cen
for(i in 1:nrow(data)){
across[i]<-0
}
for(i in 1:ncol(data)){
	if(i<6)
	{
		temp<-((data[,i]-Cen_1)^2)/96
		across<-across + temp
	}
	if(i<11 &&i>5)
	{
		temp<-((data[,i]-Cen_2)^2)/96
		across<-across + temp

	}
	if(i<16 && i>10)
	{
		temp<-((data[,i]-Cen_3)^2)/96
		across<-across + temp

	}
	if(i<21 && i>15)
	{
		temp<-((data[,i]-Cen_4)^2)/96
		across<-across + temp

	}
	if(i<26 && i>20)
	{
		temp<-((data[,i]-Cen_5)^2)/96
		across<-across + temp

	}
	if(i<31 && i>25)
	{
		temp<-((data[,i]-Cen_6)^2)/96
		across<-across + temp
	
	}
	if(i<36 && i>30)
	{
		temp<-((data[,i]-Cen_7)^2)/96
		across<-across + temp
	}
	if(i<41 &&i>35)
	{
		temp<-((data[,i]-Cen_8)^2)/96
		across<-across + temp

	}
	if(i<46 && i>40)
	{
		temp<-((data[,i]-Cen_9)^2)/96
		across<-across + temp

	}
	if(i<51 && i>45)
	{
		temp<-((data[,i]-Cen_10)^2)/96
		across<-across + temp

	}
	if(i<56 && i>50)
	{
		temp<-((data[,i]-Cen_11)^2)/96
		across<-across + temp

	}
	if(i<61 && i>55)
	{
		temp<-((data[,i]-Cen_12)^2)/96
		across<-across + temp
	
	}
	if(i<66 && i>60)
	{
		temp<-((data[,i]-Cen_13)^2)/96
		across<-across + temp
	}
	if(i<71 &&i>65)
	{
		temp<-((data[,i]-Cen_14)^2)/96
		across<-across + temp

	}
	if(i<76 && i>70)
	{
		temp<-((data[,i]-Cen_15)^2)/96
		across<-across + temp

	}
	if(i<81 && i>75)
	{
		temp<-((data[,i]-Cen_16)^2)/96
		across<-across + temp

	}
	if(i<86 && i>80)
	{
		temp<-((data[,i]-Cen_17)^2)/96
		across<-across + temp

	}
	if(i<91 && i>85)
	{
		temp<-((data[,i]-Cen_18)^2)/96
		across<-across + temp
	
	}
	if(i<96 && i>90)
	{
		temp<-((data[,i]-Cen_19)^2)/96
		across<-across + temp
	}
	if(i<101 &&i>95)
	{
		temp<-((data[,i]-Cen_20)^2)/96
		across<-across+ temp

	}
	if(i<106 && i>100)
	{
		temp<-((data[,i]-Cen_21)^2)/96
		across<-across + temp

	}
	if(i<111 && i>105)
	{
		temp<-((data[,i]-Cen_22)^2)/96
		across<-across + temp

	}
	if(i<116 && i>110)
	{
		temp<-((data[,i]-Cen_23)^2)/96
		across<-across + temp

	}
	if(i<121 && i>115)
	{
		temp<-((data[,i]-Cen_24)^2)/96
		across<-across + temp
	
	}


}
F1<-between/across


library(gplots)
h1 <- hist(F1,breaks=seq(0,22675,0.1),plot=FALSE)
h2 <- hist(F,breaks=seq(0,22675,0.1),plot=FALSE)
hh <- cbind(h1$density,h2$density)
breaks<-seq(0,22675,0.1)
breaks<-breaks[1:226750]
rownames(hh)<-breaks
colnames(hh) <- c("RMA","fRMA")
mp <-barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA F-statistics",xlim=c(0,400),ylab="Counts",xlab="F(Between group variability/across group variability )",plot.grid=TRUE)
#arrows(x0=2.29,y0=0,y1=h2$density[-1],col="red",length=0.1)
#box()

Top THREE QUARTILES
#vec<-which(Cen>4.339362)
#vecF<-which(Cen>4.577038)
#comb_vec<-intersect(vec,vecF)
#F<-F[comb_vec]
#F1<-F1[comb_vec]
#h1 <- hist(F1,breaks=seq(0,22516,0.1),plot=FALSE)
#h <- hist(F,breaks=seq(0,22516,0.1),plot=FALSE)
#h2 <- hist(F,breaks=seq(0,22516,0.1),plot=FALSE)
#hh <- cbind(h1$density,h2$density)
#colnames(hh) <- c("RMA","fRMA")
#mp <-barplot2(t(hh),beside=TRUE,col=c("red","blue"),legend=c("RMA","fRMA"),main="RMA vs fRMA F-statistics",xlim=c(0,400),ylab="Counts",plot.grid=TRUE)




